#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "glib.h"
#include "cosmoparam.h"

#define NOPERIOD 0

#define ERRFILE 1
#define ERRUSAG 2

void grid(float *x,float *y,float *z, float rlbox,  int ngrid,
	  int nobj, GHashTable * lirst, int *ll);

void busv(float *x, float *y, float *z, int ic, float rlbox, 
	  int ngrid, int nobj, GHashTable * lirst, 
	  int *ll, int *grupo, float r0, 
	  int iper, int *listvec, int *nvec);

void iden(float *x, float *y, float *z, float rlbox, float r0, 
	  int nobj, int iper, int *grupo, int *numgrupos);

void leeheader(char *nombrefile);
void lee(char *filename,float *x,float *y,float *z,int *id,int *ind);
double deltavir( double z );

struct io_header_1
{
  int      npart[6];
  double   mass[6];
  double   time;
  double   redshift;
  int      flag_sfr;
  int      flag_feedback;
  int      npartTotal[6];
  int      flag_cooling;
  int      num_files;
  double   BoxSize;
  double   Omega0;
  double   OmegaLambda;
  double   HubbleParam; 
  char     fill[256 - 6*4 - 6*8 - 2*8 - 2*4 - 6*4 - 2*4 - 4*8];  /* fills to 256 Bytes */
} header1;

int main(int argc, char **argv)
{
  int   i, l, gn, ind;
  int nobj;
  float *x,*y,*z;
  int *id, *npgrupo;
  float rlbox, redshift;
  int iper;
  int *grupo;
  double r0;
  char filename[200];
  FILE *pfout;
  int k;
  int numerodegrupos, ifile;
  int *cabecera, *lista;
  int files = 1;

  /* unused variables, should we remove them? */
  /* int dummy, n, pc; */
  /* FILE  *pf; */
  /* float r[3], v[3]; */
  /* float *vxmed, *vymed, *vzmed, *sigma; */ 

  if (argc != 3) {
    printf("usage: iden <file> <outfile>\n");
    exit(ERRUSAG);
  }

  /* READ FILE HEADER */

  if (files>1)
    sprintf(filename, "%s.0", argv[1]);
  else
    sprintf(filename, "%s", argv[1]);

  printf("Reading file: %s\n", filename);

  leeheader(filename);

  nobj = 0;
  for(k = 0; k < 6; k++)
    nobj += header1.npartTotal[k];

  rlbox=header1.BoxSize;
  redshift=header1.redshift;

  printf("Numero de particulas: %d Box en kpc: %f redshift: %f\n",
  	 nobj, header1.BoxSize, header1.redshift);
  printf("Omegas %f %f\n", header1.Omega0, header1.OmegaLambda);

  cp.omegam = header1.Omega0;
  cp.omegal = header1.OmegaLambda;
  cp.omegak = 1. - cp.omegam - cp.omegal;

  /* END OF READING FILE HEADERS */

  x          = malloc(nobj*sizeof(float));
  y          = malloc(nobj*sizeof(float));
  z          = malloc(nobj*sizeof(float));

  id         = malloc(nobj*sizeof(int));
  grupo      = malloc(nobj*sizeof(int));


  for (ifile = 0, ind = 0; ifile < files; ifile++)
  {
    printf("file ----- %d\n",ifile);
    if (files > 1)
      sprintf(filename, "%s.%d", argv[1], ifile);
    else
      sprintf(filename, "%s", argv[1]);
    printf("%s\n",filename);
    lee(filename, x, y, z, id, &ind);
  }
  
  
  iper = 1;
  r0 = rlbox/pow((float) nobj, 0.33333) *
       pow(1.0/(0.4*deltavir(redshift)), 0.33333);

  printf("nobj rlbox r0 %d %f %f\n",nobj,rlbox,r0);
  iden(x, y, z, rlbox, r0, nobj, iper, grupo, &numerodegrupos);
  printf("numero de grupos: %d\n",numerodegrupos);

  free(x);
  free(y);
  free(z);

  cabecera   = malloc(numerodegrupos*sizeof(int));
  lista      = malloc(nobj*sizeof(int));
  npgrupo    = malloc(numerodegrupos*sizeof(int));

  for (i = 0; i < numerodegrupos; i++) {
    cabecera[i] = 0;
    npgrupo[i] = 0;
  }

  for (i = 0; i < nobj; i++)
    lista[i] = 0;

  for (i = 0; i < nobj; i++) {
    cabecera[grupo[i]] = i;
    npgrupo[grupo[i]]++;
  }

  for(i=0;i<nobj;i++)
  {
    lista[cabecera[grupo[i]]]=i;
    cabecera[grupo[i]]=i;
  }

  sprintf(filename, "./data/%s.iden.dat", argv[2]);
  if ( !(pfout = fopen(filename, "w"))) {
    printf("Can't open file `%s`\n", filename);
    exit(ERRFILE);
  }


  gn=0;
  for(i=1;i<numerodegrupos;i++)
  {
    if(npgrupo[i]<10)continue;
    fwrite(&npgrupo[i],sizeof(int),1,pfout);
    l=cabecera[i];
    do
    {
      l=lista[l];
      fwrite(&l,sizeof(int),1,pfout);
    }while( l != cabecera[i] );
    fwrite(&npgrupo[i],sizeof(int),1,pfout);
    gn++;
  }
  printf("numero total de grupos com mas de 10: %d\n",gn);
  fclose(pfout);
  
  return 0;

}

void iden(float *x, float *y, float *z, float rlbox, float r0, 
	   int nobj, int iper,  int *grupo, int *numgrupos)
{
  /*se asume que las posisiones estan entre 0 y rlbox*/
  int   i, j, k;
  GHashTable * lirst; 
  int   *ll;
  int   *listvec, *listvecnew;
  int   ngrid;
  int   nvec, nvecvec;
  int   ngrupos;

  ngrid = (int) (rlbox/r0);
  printf("BEFORE r0=%f ngrid=%d\n", r0, ngrid);
  if(ngrid > 512)
    ngrid = 512;
  printf("AFTER r0=%f ngrid=%d\n", r0, ngrid);

  /* This is a lot of memory! If nobj = 128**3 then this is
     128 ** 3 * 32 * 3 + 512 * 512 * 512 * 32 ~ 536MB ! */
  lirst      = g_hash_table_new_full(g_int_hash, g_int_equal, free, free);

  ll         = malloc(nobj*sizeof(int));
  listvec    = malloc(nobj*sizeof(int));
  listvecnew = malloc(nobj*sizeof(int));

  /* Construct the grid. */
  grid(x, y, z, rlbox, ngrid, nobj, lirst, ll);

  /* Empieza la identificacion */

  /* Initialization */
  for( i = 0 ; i < nobj ; i++ )
    grupo[i] = 0;

  printf("Empieza Iden\n");
  fflush(stdout); 

  ngrupos = 0;
  for( i = 0 ; i < nobj; i++ )
  {
    if (i % 100000 == 0) {
      printf("%i\r",i); 
      fflush(stdout);
    }
    
    if (grupo[i] != 0) 
      continue; /* salta a la siguiente */
    nvec = 0;
    busv(x, y, z, i, rlbox, ngrid, nobj, lirst, ll, grupo, r0, iper, listvec, 
	 &nvec);

    if(nvec!=0)
    {
      ngrupos++;
      grupo[i] = ngrupos;
    }
    for( j = 0 ; j < nvec ; j++ )
      grupo[listvec[j]] = ngrupos;

    do
    {
      nvecvec=0;
      for (j = 0; j < nvec; j++)
      {
        busv(x, y, z, listvec[j], rlbox, ngrid, nobj, lirst, ll, grupo, r0, iper,
	     listvecnew, &nvecvec);
        for (k = 0; k < nvecvec; k++)
          grupo[listvecnew[k]] = ngrupos;
      }

      if( nvecvec == 0 ) break;
      for (j = 0 ; j < nvecvec ; j++)
        listvec[j]=listvecnew[j];
      nvec=nvecvec;
    } while( 1 ); /*fin del do de los vecinos de los vecinos*/
  } 

  *numgrupos = ngrupos;
  printf("termino identificacion\n");

  /* termino la identificacion */

  free(ll);
  g_hash_table_destroy(lirst);
  free(listvec);
  free(listvecnew);
}

void grid(float *x,float *y,float *z, float rlbox, int ngrid, 
	  int nobj, GHashTable * lirst, int *ll)
{

  int i, j, k;
  int ix, iy, iz;
  float fac;

  /* Each box in the grid (each index in lrist) represents
     "ngrid/rlbox" volume of space. That is each box is "ngrid/rlbox"
     long by "ngrid/rlbox" width by "ngrid/rlbox" high. */
   fac = (float) ngrid/rlbox ;

  /* Initialize list of objects. */
  for (i = 0 ; i < nobj ; i++ )  
    ll[i] = 0 ;

  /* Fill the grid (lirst) with one object each. Several objects fit
     in one box of the grid, that is several objects have the same ix,
     iy, iz coordinates (when we cast the x, y, z coordinates to
     ints). The object which we keep in a box with coordinates ix, iy,
     iz is the last object from the list (enumarting the objcets from
     0 to nobj) that has those coordinates (after the cast). */
  int * key, * val;
  for (i = 0 ; i < nobj ; i++ ) {
    ix = (int) (x[i]*fac);
    iy = (int) (y[i]*fac);
    iz = (int) (z[i]*fac);
    key = malloc(sizeof(int));
    *key = (ix * ngrid + iy) * ngrid + iz;
    val = malloc(sizeof(int));
    *val = i;
    g_hash_table_insert(lirst, key, val);
  }  

  /* With this loop we create, for every box with ix, iy, iz,
     coordinates, a list of the objects that should be in a box.
 
     Say that x_j is the index of the last object that was placed in a
     box of the grid (lirst) by the loop above. This loop creates the
     list in ll by having the x_j index of ll pointing to the first
     object (from 0 to nobj) which should be in ix, iy, iz box and
     this object points to the second and so forth such that 
     x_(j-1) points back to x_j again.

     ll   :    |x1| | |x2| | |     | |  x_j  |  | |  | x0|
     index:    |x0| | |x1| | | ... | |x_(j-1)|  | |  |x_j|
     
     All list are place in ll.
  */

  gpointer * idx;
  for( i = 0 ; i < nobj ; i++ )
  {
    ix = (int) (x[i]*fac);
    iy = (int) (y[i]*fac);
    iz = (int) (z[i]*fac);
    key = malloc(sizeof(int));
    *key = (ix * ngrid + iy) * ngrid + iz;
    idx = g_hash_table_lookup(lirst, key);
    ll[(idx != NULL ? (int) * idx : 0)] = i;
    val = malloc(sizeof(int));
    *val = i;
    g_hash_table_replace(lirst, (gpointer) key, (gpointer) val);
  }   
  
  return;
}

void busv(float *x, float *y, float *z, int ic, float rlbox, 
	  int ngrid, int nobj, GHashTable *lirst, 
	  int *ll, int *grupo, float r0, int iper, 
	  int *listvec, int *nvec)
{
  int ixc, iyc, izc;
  int ixci, iyci, izci;
  int ixcf, iycf, izcf;
  int ix, iy, iz;
  int ixx, iyy, izz;
  int ibox;
  float xx, yy, zz, dis;
  int i;
  float rlbox2;

  rlbox2=rlbox/2.0;

  ixc = (int) (x[ic]/rlbox*ngrid);
  ixci = ixc-1;
  ixcf = ixc+1;

  iyc=(int)(y[ic]/rlbox*ngrid);
  iyci=iyc-1;
  iycf=iyc+1;

  izc=(int)(z[ic]/rlbox*ngrid);
  izci=izc-1;
  izcf=izc+1;

  /* ixc, iyc, izc is the coordinate of the box for the object given
     by the coordinates x[ic], y[ic], z[ic]. We also need the
     coordinates of the boxes next to this one. The coordinates of the
     nieghbor boxes are given by conbinations of ixci, ixcf, iyci,
     iycf, izci, izcf. That's what are the "for loops" for, iterating
     over these boxes. */
  
  if ( iper == NOPERIOD )
  {
    if( ixci < 0 ) ixci = 0;
    if( iyci < 0 ) iyci = 0;
    if( izci < 0 ) izci = 0;
    if( ixcf >= ngrid ) ixcf = ngrid-1;
    if( iycf >= ngrid ) iycf = ngrid-1;
    if( izcf >= ngrid ) izcf = ngrid-1;
  } 

  for (ixx = ixci ; ixx <= ixcf ; ixx++)
  {
    ix = ixx;

    if (iper == 1) {
      if (ix >= ngrid)
	ix = ix - ngrid;
      
      if (ix < 0) 
	ix = ix + ngrid;
    }

    for (iyy = iyci ; iyy <= iycf ; iyy++) {
      iy = iyy;
      
      if (iper == 1) {
        if(iy >= ngrid)
	  iy = iy - ngrid;
        if(iy < 0)
	  iy = iy + ngrid;
      }
      
      for( izz = izci ; izz <= izcf ; izz++) {
        iz = izz;
        if(iper == 1) {
          if(iz >= ngrid)
	    iz = iz - ngrid;
          if(iz < 0)
	    iz = iz + ngrid;
        }


        ibox = (ix * ngrid + iy) * ngrid + iz ;
       
        gpointer * last;
	int last_value;
        last = g_hash_table_lookup(lirst, (gpointer) &ibox);
	last_value = last != NULL ? (int) * last : 0; 
	i = last_value;

        if (i == -1)
	  continue;

	do
        {
	  /* Get the first object of the list. Note that lirst[index]
	     has the "last" object which belongs to that box and the
	     "last" object points to the "first" in ll */
          i = ll[i];

	  /* Calculate the distance between the object with index "ic"
	     and the next object given by index "i". */
          xx = x[i]-x[ic];
          yy = y[i]-y[ic];
          zz = z[i]-z[ic];
          
	  /* Take into account perdiodicy. */
	  if (iper)
          {
            if( xx > rlbox2 ) xx = xx - rlbox;
            if( yy > rlbox2 ) yy = yy - rlbox;
            if( zz > rlbox2 ) zz = zz - rlbox;
            if( xx < -rlbox2 ) xx = xx + rlbox;
            if( yy < -rlbox2 ) yy = yy + rlbox;
            if( zz < -rlbox2 ) zz = zz + rlbox;
          }
          
	  dis=sqrt(xx*xx+yy*yy+zz*zz);
	  
	  /* if the distance is within "friend distance" add it to the
	     friends list */
          if( dis < r0 && grupo[i] == 0 && i != ic )
          {
            listvec[*nvec]=i;
            *nvec=*nvec+1;
          }
	  /* last = g_hash_table_lookup(lirst, (gpointer) &ibox); */
	  /* value = (last != NULL ? (int) * last : 0); */
	  /* printf("%d, %d\n", i, value); */
	  /* next = g_hash_table_lookup(lirst, (gpointer) &ibox); */
        } while (i != last_value); /* cycle through the circular
				      list formed by ll */

      } /*fin izz*/
  
    } /*fin iyy*/
  
  } /*fin ixx*/

}

void leeheader(char *filename)
{
  FILE *pf;
  int dummy;

  if ( !(pf = fopen(filename, "r"))) {
    printf("Can't open file `%s`\n", filename);
    exit(ERRFILE);
  }

  fread(&dummy, sizeof(dummy), 1, pf);
  fread(&header1, sizeof(header1), 1, pf);
  fread(&dummy, sizeof(dummy), 1, pf);
  fclose(pf);
}

void lee(char *filename, float *x, float *y, float *z, int *id, int *ind)
{
  FILE *pf;
  int dummy, k, pc, n;
  float r[3],v[3];
  
  if ( !(pf = fopen(filename, "r"))) {
    printf("Can't open file `%s`\n", filename);
    exit(ERRFILE);
  }

  fread(&dummy, sizeof(dummy), 1, pf);
  fread(&header1, sizeof(header1), 1, pf);
  fread(&dummy, sizeof(dummy), 1, pf);
  fread(&dummy, sizeof(dummy), 1, pf);
  for(k=0,pc=0;k<6;k++)
  {
    for(n=0;n<header1.npart[k];n++)
    {
      fread(&r[0], sizeof(float), 3, pf);
      if(r[0] < 0.)    r[0] += header1.BoxSize ;
      if(r[0] >= header1.BoxSize ) r[0] -= header1.BoxSize  ;
      if(r[1] < 0.)    r[1] += header1.BoxSize  ;
      if(r[1] >= header1.BoxSize ) r[1] -= header1.BoxSize  ;
      if(r[2] < 0.)    r[2] += header1.BoxSize  ;
      if(r[2] >= header1.BoxSize ) r[2] -= header1.BoxSize  ;
      x[*ind+pc]=r[0];
      y[*ind+pc]=r[1];
      z[*ind+pc]=r[2];

      if(z[pc]<1.0e-10) 
        printf("------------------ %g %g %g\n",x[pc],y[pc],z[pc]);
      pc++;
    }
  }
  fread(&dummy, sizeof(dummy), 1, pf);
  printf("pc pos %d\n",pc);
  fread(&dummy, sizeof(dummy), 1, pf);
  for(k=0;k<6;k++)
  {
    for(n=0;n<header1.npart[k];n++)
    {
      fread(&v[0], sizeof(float), 3, pf);
      pc++;
    }
  }
  fread(&dummy, sizeof(dummy), 1, pf);
  fread(&dummy, sizeof(dummy), 1, pf);
  for(k=0,pc=0;k<6;k++)
  {
    for(n=0;n<header1.npart[k];n++)
    {
      fread(&id[*ind+pc], sizeof(int), 1, pf);
      pc++;
    }
  }
  printf("pc ind %d\n",pc);

  *ind+=pc;
  printf("--------------------- ind : %d\n",*ind);
  fclose(pf);
}

double deltavir( double z )
{
    double q ;
    double as, x ;
    as=1./(1.+z);
    if ((cp.omegam == 1.) && (cp.omegal == 0.))
    {
        q = 177.7 ;      
        /* Einstein - de Sitter Model Peebles, 1980. Eq 19.50 */
    }
    if ((cp.omegam < 1.) && (cp.omegal == 0.))
    {
        x = as * 2. *(1./cp.omegam-1.);
        q = pow(2.*M_PI,2.)*pow(x,3.)/pow(sqrt(x*(x+2.))-log((1.+x)+sqrt(x*(x+2.))),2.);
        /* Maoz, 1990 Eq 12. Lacey & Cole, 1993 Eq A15. Oukbir & Blanchard, 1997 */
    }
    if ((cp.omegam < 1.) && (cp.omegak == 0.))
    {
        x = as*pow(1./cp.omegam-1.,1./3.);
        q = 18.*M_PI*M_PI*(1.+0.4093*pow(x,2.71572)); 
        /* Nakamura & Suto, 1997 Eq c-19. Eke et al, 1996 */
    }
    return q;
}
